home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / linpkdrv.zip / SP.FOR < prev    next >
Text File  |  1984-01-07  |  17KB  |  613 lines
  1. C     MAIN PROGRAM
  2.       INTEGER LUNIT
  3. C     ALLOW 5000 UNDERFLOWS.
  4. C     CALL TRAPS(0,0,5001,0,0)
  5. C
  6. C     OUTPUT UNIT NUMBER
  7. C
  8.       LUNIT = 6
  9. C
  10.       CALL SPOTS(LUNIT)
  11.       STOP
  12.       END
  13.       SUBROUTINE SPOTS(LUNIT)
  14. C     LUNIT IS THE OUTPUT UNIT NUMBER
  15. C
  16. C     TESTS
  17. C        SPOCO,SPOFA,SPOSL,SPODI,SPPCO,SPPFA,SPPSL,SPPDI
  18. C        SPBCO,SPBFA,SPBSL,SPBDI
  19. C
  20. C     LINPACK.  THIS VERSION DATED 08/14/78 .
  21. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  22. C
  23. C     SUBROUTINES AND FUNCTIONS
  24. C
  25. C     LINPACK SPOCO,SPOSL,SPODI,SPPCO,SPPSL,SPPDI
  26. C     LINPACK SPBCO,SPBSL,SPBDI
  27. C     EXTERNAL SPOXX,SMACH
  28. C     BLAS SAXPY,SDOT,SASUM
  29. C     FORTRAN ABS,AMAX1,FLOAT,MAX0
  30. C
  31. C     INTERNAL VARIABLES
  32. C
  33.       REAL A(15,15),AB(15,15),AINV(15,15),ASAVE(15,15)
  34.       REAL AP(120),B(15),SDOT,X(15),XB(15),XEXACT(15)
  35.       REAL XP(15),T,Z(15)
  36.       REAL ANORM,AINORM,COND,COND1,DET(2),DETB(2),DETP(2)
  37.       REAL EN,ENORM,EPS,FNORM,ONEPX,Q(6),QS(6),RCOND,RCONDB
  38.       REAL RCONDP,RNORM,S,SASUM,SMACH,XNORM
  39.       INTEGER I,INFO,INFOB,INFOP,IQ(6),I1,J,JB
  40.       INTEGER K,KASE,KB,KBFAIL,KNPD,KOUNT,KPFAIL
  41.       INTEGER KSUSP(6),LDA,LUNIT,M,N,NPRINT
  42.       LOGICAL KBF,KPF
  43. C
  44.       LDA = 15
  45. C
  46. C     WRITE MATRIX AND SOLUTIONS IF  N .LE. NPRINT
  47. C
  48.       NPRINT = 3
  49. C
  50.       WRITE (LUNIT,560)
  51.       WRITE (LUNIT,1000)
  52. C
  53.       DO 10 I = 1, 6
  54.          KSUSP(I) = 0
  55.    10 CONTINUE
  56.       KNPD = 0
  57.       KPFAIL = 0
  58.       KBFAIL = 0
  59. C
  60. C     SET EPS TO ROUNDING UNIT FOR REAL ARITHMETIC
  61. C
  62.       EPS = SMACH(1)
  63.       WRITE (LUNIT,570) EPS
  64.       WRITE (LUNIT,550)
  65. C
  66. C     START MAIN LOOP
  67. C
  68.       KASE = 1
  69.    20 CONTINUE
  70. C
  71. C        GENERATE TEST MATRIX
  72. C
  73.          CALL SPOXX(A,LDA,N,KASE,LUNIT)
  74. C
  75. C        N = 0 SIGNALS NO MORE TEST MATRICES
  76. C
  77. C     ...EXIT
  78.          IF (N .LE. 0) GO TO 540
  79.          ANORM = 0.0E0
  80.          DO 30 J = 1, N
  81.             ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
  82.    30    CONTINUE
  83.          WRITE (LUNIT,720) ANORM
  84. C
  85.          IF (N .GT. NPRINT) GO TO 50
  86.             WRITE (LUNIT,550)
  87.             DO 40 I = 1, N
  88.                WRITE (LUNIT,760) (A(I,J), J = 1, N)
  89.    40       CONTINUE
  90.             WRITE (LUNIT,550)
  91.    50    CONTINUE
  92. C
  93. C        GENERATE EXACT SOLUTION
  94. C
  95.          XEXACT(1) = 1.0E0
  96.          IF (N .GE. 2) XEXACT(2) = 0.0E0
  97.          IF (N .LE. 2) GO TO 70
  98.             DO 60 I = 3, N
  99.                XEXACT(I) = -XEXACT(I-2)
  100.    60       CONTINUE
  101.    70    CONTINUE
  102. C
  103. C        SAVE MATRIX AND GENERATE R.H.S.
  104. C
  105.          DO 90 I = 1, N
  106.             B(I) = 0.0E0
  107.             DO 80 J = 1, N
  108.                ASAVE(I,J) = A(I,J)
  109.                B(I) = B(I) + A(I,J)*XEXACT(J)
  110.    80       CONTINUE
  111.             X(I) = B(I)
  112.             XP(I) = X(I)
  113.             XB(I) = X(I)
  114.    90    CONTINUE
  115. C
  116. C        FACTOR AND ESTIMATE CONDITION
  117. C
  118.          RCOND = -1.0E0
  119.          CALL SPOCO(A,LDA,N,RCOND,Z,INFO)
  120. C
  121. C        OUTPUT NULL VECTOR
  122. C
  123.          IF (N .GT. NPRINT .OR. INFO .NE. 0) GO TO 110
  124.             WRITE (LUNIT,770)
  125.             DO 100 I = 1, N
  126.                WRITE (LUNIT,780) Z(I)
  127.   100       CONTINUE
  128.             WRITE (LUNIT,550)
  129.   110    CONTINUE
  130. C
  131. C        FACTOR PACKED FORM AND COMPARE
  132. C
  133.          KPF = .FALSE.
  134.          K = 0
  135.          DO 130 J = 1, N
  136.             DO 120 I = 1, J
  137.                K = K + 1
  138.                AP(K) = ASAVE(I,J)
  139.   120       CONTINUE
  140.   130    CONTINUE
  141.          RCONDP = -1.0E0
  142.          CALL SPPCO(AP,N,RCONDP,Z,INFOP)
  143.          IF (INFOP .EQ. INFO) GO TO 140
  144.             WRITE (LUNIT,880)
  145.             WRITE (LUNIT,920) INFO,INFOP
  146.             KPF = .TRUE.
  147.   140    CONTINUE
  148.          IF (RCONDP .EQ. RCOND) GO TO 150
  149.             WRITE (LUNIT,880)
  150.             WRITE (LUNIT,930) RCOND,RCONDP
  151.             KPF = .TRUE.
  152.   150    CONTINUE
  153.          K = 0
  154.          KOUNT = 0
  155.          DO 170 J = 1, N
  156.             DO 160 I = 1, J
  157.                K = K + 1
  158.                IF (AP(K) .NE. A(I,J)) KOUNT = KOUNT + 1
  159.   160       CONTINUE
  160.   170    CONTINUE
  161.          IF (KOUNT .EQ. 0) GO TO 180
  162.             WRITE (LUNIT,880)
  163.             WRITE (LUNIT,940) KOUNT
  164.             KPF = .TRUE.
  165.   180    CONTINUE
  166. C
  167. C        FACTOR BAND FORM AND COMPARE
  168. C
  169.          KBF = .FALSE.
  170.          M = 0
  171.          DO 200 J = 1, N
  172.             DO 190 I = 1, J
  173.                IF (ASAVE(I,J) .NE. 0.0E0) M = MAX0(M,J-I)
  174.   190       CONTINUE
  175.   200    CONTINUE
  176. C
  177.          DO 220 J = 1, N
  178.             I1 = MAX0(1,J-M)
  179.             DO 210 I = I1, J
  180.                K = I - J + M + 1
  181.                AB(K,J) = ASAVE(I,J)
  182.   210       CONTINUE
  183.   220    CONTINUE
  184.          WRITE (LUNIT,840) M
  185.          RCONDB = -1.0E0
  186.          CALL SPBCO(AB,LDA,N,M,RCONDB,Z,INFOB)
  187.          IF (INFOB .EQ. INFO) GO TO 230
  188.             WRITE (LUNIT,830)
  189.             WRITE (LUNIT,920) INFO,INFOB
  190.             KBF = .TRUE.
  191.   230    CONTINUE
  192.          IF (RCONDB .EQ. RCOND) GO TO 240
  193.             WRITE (LUNIT,830)
  194.             WRITE (LUNIT,930) RCOND,RCONDB
  195.             KBF = .TRUE.
  196.   240    CONTINUE
  197.          KOUNT = 0
  198.          DO 250 J = 1, N
  199.             IF (AB(M+1,J) .NE. A(J,J)) KOUNT = KOUNT + 1
  200.   250    CONTINUE
  201.          IF (KOUNT .EQ. 0) GO TO 260
  202.             WRITE (LUNIT,830)
  203.             WRITE (LUNIT,940) KOUNT
  204.             KBF = .TRUE.
  205.   260    CONTINUE
  206. C
  207. C        TEST FOR DEFINITENESS
  208. C
  209.          IF (INFO .EQ. 0) GO TO 270
  210.             WRITE (LUNIT,860) INFO
  211.             KNPD = KNPD + 1
  212.          GO TO 530
  213.   270    CONTINUE
  214.             COND = 1.0E0/RCOND
  215.             WRITE (LUNIT,590) COND
  216.             ONEPX = 1.0E0 + RCOND
  217.             IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,580)
  218. C
  219. C           COMPUTE INVERSE, DETERMINANT AND COND1 = TRUE CONDITION
  220. C
  221.             DO 290 J = 1, N
  222.                DO 280 I = 1, J
  223.                   AINV(I,J) = A(I,J)
  224.   280          CONTINUE
  225.   290       CONTINUE
  226.             CALL SPODI(AINV,LDA,N,DET,11)
  227.             AINORM = 0.0E0
  228.             DO 310 J = 1, N
  229.                DO 300 I = J, N
  230.                   AINV(I,J) = AINV(J,I)
  231.   300          CONTINUE
  232.                AINORM = AMAX1(AINORM,SASUM(N,AINV(1,J),1))
  233.   310       CONTINUE
  234.             COND1 = ANORM*AINORM
  235.             WRITE (LUNIT,600) COND1
  236.             WRITE (LUNIT,800) DET(1)
  237.             WRITE (LUNIT,810) DET(2)
  238. C
  239. C           SOLVE A*X = B
  240. C
  241.             CALL SPOSL(A,LDA,N,X)
  242. C
  243.             IF (N .GT. NPRINT) GO TO 330
  244.                WRITE (LUNIT,610)
  245.                DO 320 I = 1, N
  246.                   WRITE (LUNIT,790) X(I)
  247.   320          CONTINUE
  248.                WRITE (LUNIT,550)
  249.   330       CONTINUE
  250. C
  251. C           MORE PACKED COMPARE
  252. C
  253.             CALL SPPSL(AP,N,XP)
  254.             KOUNT = 0
  255.             DO 340 I = 1, N
  256.                IF (XP(I) .NE. X(I)) KOUNT = KOUNT + 1
  257.   340       CONTINUE
  258.             IF (KOUNT .EQ. 0) GO TO 350
  259.                WRITE (LUNIT,880)
  260.                WRITE (LUNIT,950) KOUNT
  261.                KPF = .TRUE.
  262.   350       CONTINUE
  263.             CALL SPPDI(AP,N,DETP,11)
  264.             IF (DETP(1) .EQ. DET(1) .AND. DETP(2) .EQ. DET(2))
  265.      *         GO TO 360
  266.                WRITE (LUNIT,880)
  267.                WRITE (LUNIT,960) DETP
  268.                KPF = .TRUE.
  269.   360       CONTINUE
  270.             KOUNT = 0
  271.             K = 0
  272.             DO 380 J = 1, N
  273.                DO 370 I = 1, J
  274.                   K = K + 1
  275.                   IF (AP(K) .NE. AINV(I,J)) KOUNT = KOUNT + 1
  276.   370          CONTINUE
  277.   380       CONTINUE
  278.             IF (KOUNT .EQ. 0) GO TO 390
  279.                WRITE (LUNIT,880)
  280.                WRITE (LUNIT,970) KOUNT
  281.                KPF = .TRUE.
  282.   390       CONTINUE
  283. C
  284. C           MORE BAND COMPARE
  285. C
  286.             CALL SPBSL(AB,LDA,N,M,XB)
  287.             KOUNT = 0
  288.             DO 400 I = 1, N
  289.                IF (XB(I) .NE. X(I)) KOUNT = KOUNT + 1
  290.   400       CONTINUE
  291.             IF (KOUNT .EQ. 0) GO TO 410
  292.                WRITE (LUNIT,830)
  293.                WRITE (LUNIT,950) KOUNT
  294.                KBF = .TRUE.
  295.   410       CONTINUE
  296.             CALL SPBDI(AB,LDA,N,M,DETB)
  297.             IF (DETB(1) .EQ. DET(1) .AND. DETB(2) .EQ. DET(2))
  298.      *         GO TO 420
  299.                WRITE (LUNIT,830)
  300.